A multi-state dynamic process confers mechano-adaptation to a biological nanomachine

Adaptation is a defining feature of living systems. The bacterial flagellar motor adapts to changes in the external mechanical load by adding or removing torque-generating (stator) units. But the molecular mechanism behind this mechano-adaptation remains unclear. Here, we combine single motor eletrorotation experiments and theoretical modeling to show that mechano-adaptation of the flagellar motor is enabled by multiple mechanosensitive internal states. Dwell time statistics from experiments suggest the existence of at least two bound states with a high and a low unbinding rate, respectively. A first-passage-time analysis of a four-state model quantitatively explains the experimental data and determines the transition rates among all four states. The torque generated by bound stator units controls their effective unbinding rate by modulating the transition between the bound states, possibly via a catch bond mechanism. Similar force-mediated feedback enabled by multiple internal states may apply to adaptation in other macromolecular complexes.


A plausible model for stator dynamics and its simplified form
Here we present a conceptual model to described the stochastic motion of a bound stator unit to demonstrate a plausible mechanism for the torque-dependence of the off rate. The basic idea is that a bound (tethered) stator unit can become unbound with a unbinding rate k off (x) that can depend on its displacement x. This dependence of k off on x is reasonable as there are barriers between the bound and unbound states that can depend on location of the stator unit relative to the landing point (the attachment point at x = 0). We assume that k off (x) is maximum at x = 0 as the stator unit can unbind through the same pathway as the binding process at x = 0 and it decreases as |x| increases, e.g., k off = k 0 exp(−|x|/x 0 ) or k off = k 0 exp −x 2 /x 2 0 . This symmetric form of k off for x → −x is inspired from the CCW and the CW symmetry in stator remodeling 1 .
For a tethered stator unit, we can write down the dynamics of its displacement x as a Langevin equation: where ξ is the viscosity; η is the thermal noise with η = 0 and η(t)η(t ) = 2k B T ξ −1 δ (t − t ) (k B T is the thermal energy unit); V 0 (x) is the free energy of the tethered stator. If we assume the tethers can be approximated as springs with a spring constant k, we have: V 0 (x) = 1 2 kx 2 . Of course, the free energy V 0 (x) can have more complex landscape with multiple minima separated by energy barriers.
The driving force F = Γ/R is due to the counter force from the rotor to the stator unit with Γ the torque and R the rotor radius. Over a time scale much longer than the time scale of proton translocation through the stator units, Γ may be considered constant. We can then solve Eq. 1 and obtain the steady state distribution of x as: where A is normalization constant. Given the steady state distribution of x given in Eq. 2 and the dependence of the off rate on x given by k off (x) to determine the observed overall off ratek off : which depends on Γ because P s (x) depends on Γ.
For simplicity, we can coarse grain the continuous x-dependence by two states, T and L, where the two states are separated by a displacement threshold x th . L represents a loose state with |x| ≤ x th around x 0 where k off (x) ≡ k off,l is large, and T represents a tight state at |x| > x th where k off (x) ≡ k off,t is much smaller than k off,l . This simplified model of the bound stator is used in this work.

Dwell time statistics in two-state model with only one bound state
If there is only one bound state with a single off rate k off and an on rate k on , the model becomes much simpler. We call this simplified model the two-state model. The survival probability for state-N can be obtained analytically: where k + (N) = (N tot − N)k on is the total rate of increasing N by one with k on being the on rate for each binding site, and k − (N) = Nk off is the total rate of decreasing N by one with k off being the off rate for each bound stator.
From the distribution of survival probability, we can obtain the distribution of dwell times which follows an exponential decay function with a time scale τ ts (N) for the two-state model: τ ts depends on N since both k + (N) = (N tot − N)k on and k − (N) = Nk off depend on N. However, the range of change in τ ts (N) is limited. The ratio between the time scales for any two levels N 1 and N 2 (< N 1 ) is bounded by:

2/12
where the equality only exists in the extreme case when k off = 0. For example, with N tot = 12, the change in time scale in the two-state model between N 1 = 7 and N 2 = 2 is less than 2: τ ts (7) τ ts (2) ≤ 10 5 = 2. Fig. 3 shows the experimentally determined distribution of dwell times P exp (τ|N) for different N. The average dwell time τ (N) for a given N is then computed from P exp (τ|N). As shown in Fig. 3a in the main text, the experimentally measured timescales τ (N) change by over 3-fold as N increases from its low values (N = 1, 2, 3) to high values (N = 6, 7, 8). For example, τ (7) = 43.2 s and τ (2) = 6.6 s, so τ (7) τ (2) ≈ 6.5, which is much larger than the ratio expected from a two-state model. The disagreement between the two-state model and experiments can also be seen in the normalized variance magnitude Using Eqn. 4, one can show that V = 1 for a two-state model, whereas V ≈ 2 was observed in the experiments (Fig. 3c in the main text). Furthermore, with only one bound state, the overall off rate k − (N) = Nk off,l becomes independent of time. The distribution functions for the dwell time followed by adding ('+') or removing ('-') of a stator can then be determined analytically for the two-state model: In this case, the average dwell times for the '+' (on) and the '-' (off) transitions are equal: τ + (N) = τ − (N) = (k − (N) + k + (N)) −1 , which is again inconsistent with experiments, which show that τ + (N) = τ − (N) (Fig. 4).
Taken together, our experimental results rule out the simple two-state model with only one bound state and suggest the existence of multiple bound states with multiple time scales in the remodeling process.

Dwell time statistics with two bound states
Here, we determine the dwell time statistics in our proposed model ( Fig. 1b in the main text) with two bound states (L and T) by using first-passage time analysis. In this analysis, we neglect the H-state that only affects the short time behavior of the system.
Assume at t = 0, the motor enters into the state-N with N engaged stator units (by either adding or removing a stator unit).
The motor exits the state-N by either having a new stator unit engaged with the rotor (the "on" process) or by having one of the N bound stator units leaving the motor (the "off" process), whichever happens first.
Since the observed transition (either on or off) depends on which transition happens first, the problem is a first passage time problem. The first passage time problem can be treated by computing the survival probability S(t), which is the probability the motor stays in state-N at time t ≥ 0. Obviously, we have S(0) = 1 and S(∞) = 0.

3/12
There are two processes (on and off) that terminate the state-N. The on process has a constant (time-independent) rate: where N tot is the maximum number of stator units each motor can accommodate, and the on rate k on can depend on the rotational speed ω, which is linearly proportional to the number of stator units (ω = ω 1 N in the high load regime with ω 1 the rotational speed per stator unit). Additionally, k on can also depend directly on N due to other mechanisms, such as co-operative interactions.
For N ≥ 1, the off process has a rate k − (t) that is more complex as there are multiple stator units each with a different effective "off" rate depending on their arrival time t N < t N−1 < .. < t 1 ≤ 0 where the newest stator unit arrives at t 1 ≤ 0. For a given unit-i, the probability of it being still bound with the rotor is its "survival" probability S i (t). The survived unit-i can be either in the L or the T state with probability L i (t) and T i (t), respectively, and S i = T i + L i . L i (t) and T i (t) can be determined by solving the following equations: with the initial conditions L i (t i ) = 1 and T i (t i ) = 0 as we assume the stator is in the L state when it initially binds to the rotor.
Since a stator unit has a much smaller off rate in the T state as compared to the L state, i.e., k off,t << k off,l , by setting k off,t = 0 we can approximate the overall off rate k − as: where p i (t) = L i L i +T i is the fractional probability of unit-i to be in the L state at time t.
Given the rates k + and k − , the survival probability satisfies: S(t + dt) = S(t) × (1 − (k + + k − )dt). By taking the limit 4/12 dt → 0, we have: which can be solved with the initial condition S(0) = 1: which is the product of the survival probabilities of all the individual bound stator units and the survival probability (exp(−k + t)) due to adding a stator unit. Note that in the product we have S i (t) S i (0) to satisfy S(0) = 1.
From Eq. 12, the expression for S i can be written as: where c = σ + −k l −k t σ + −σ − is a constant. From the expression of S i , we can obtain the expression for S i (t) S i (0) : where the coefficient c i is given by: where ∆σ ≡ σ + − σ − > 0. It is easy to see that c i < 1 and since c can be negative, c i can also be negative. Given the order of the stator unit arrival times: t N < t N−1 < ... < t 1 ≤ 0, the absolute value |c i | decreases with i (in fact, |c i | decreases exponentially with |t i |).
Finally, the expression for S(t|N) can be written as: which contains (N + 1) exponential decay terms each with a different decay rate (or timescale). The coefficients are: Given that c i can be negative, a 1 can be negative. In principle, a n can depend on N.
For simplicity, we can keep only the larger values of c i with i ≤ i m and let the other values of c i to be zero c i>i m = 0. For

5/12
example, the simplest model corresponds to i m = 1 and we have: where k a = k + + Nσ − , k b = k a + ∆σ , and S(t) = e −k +,0 t for N = 0. The model then contains three parameters σ ± and c 1 ≈ c if we assume t 1 ≈ 0.
For the results reported in Figure 3 in the main text, we used i m = 3, i.e., by making the approximation that c i≥4 = 0 and treating c 1 , c 2 , and c 3 as fitting parameters. This leads to an explicit expression for S(t|N): where k c = k a + 2∆σ and k d = k a + 3∆σ . The normalization condition S(0) = 1 is preserved in the expression for S(t) given above, which has four decay rates (terms) when N ≥ 3 and there are three variables c 1 , c 2 , and c 3 that we used to fit the experimental data. Given that the arrival time for the last stator unit t 1 = 0 (at least for most cases), we have: c 1 ≈ c. The value of c reported in the main text is the same as c 1 determined by fitting the model to the experiments.
From the expression of S(t|N) given in Eq. 21, we can obtain the expressions for k − (τ|N) = − d ln S dt | t=τ − k + and the two dwell time distributions P + (τ|N) and P − (τ|N) for the '+' and '-' transitions when the bound stator number is N: The overall dwell-time distribution P(τ|N) = P + (τ|N) + P − (τ|N). These dwell time distributions can be used to derive expressions for τ , f + , and V for each N: which depend on k + (N), σ − , ∆σ , c 1 , c 2 and c 3 .

Possible role of cell-to-cell variability
Here we consider the possible effects of cell-to-cell variability on the statistics of dwell times and in particular on the variability in dwell times for a given stator unit number N. If there is a large enough variability in rates, especially k + due to cell-to-cell variation in the number of freely diffusing stator units, then even within a 2-state model such cell-to-cell variability could potentially account for the large observed values of V ≡ τ 2 τ 2 − 1, where τ is the dwell time. Such cell-to-cell variability would increase V . However, as we show below, the large value of V (≈ 2) observed in our experiments can not be explained by cell-to-cell variability because that would require an unrealistically high level of variability in the kinetic rates in a highly correlated fashion.
In a 2-state model, the statistics of dwell time are determined by the total rate k = k + + k − , which is the sum of the on-rate (k + ) and the off-rate (k − ). For a given k, τ follows an exponential distribution: P(τ|k) = k exp(−kτ) with its first and second moments given by τ = k −1 and τ 2 = 2k −2 , respectively. Now, let us assume that the cell-to-cell variation in k can be described by a distribution q(k) for k. We can then express the normalized variance V as: For simplicity, let us assume that k has a broad uniform distribution between k 1 and k 2 , i.e., q(k) = (k 2 − k 1 ) −1 when k 1 ≤ k ≤ k 2 and 0 otherwise. This distribution is characterized by a maximum fold change f = k 2 /k 1 ≥ 1. From Equation 26, one can show that V depends on the fold change f as: which has a rather weak dependence on f (see Fig. 5 above). When f → 1, i.e., there is no variation in k, V → 1. V increases with f , but slowly. For V to reach the observed value of ≈ 2, f has to be as large as 10. For a more reasonable value of f = 4, we have V = 1.34, which is far below the observed value of V .
Another important issue is that since k = k + + k − , a large fold change in k requires both k + and k − to have the same (or similar) fold change in a correlated way. If the variations in k + and k − are uncorrelated, which we believe is the case here, the variation in k measured by fold change will be dominated by the rate that has the higher average value. At higher value of N, k − is comparable to and even larger than k + , in which case the variation in k could be dominated by variation in k − . However, within the 2-state model, there is no obvious reason for k − to have a large variability. For example, if we assume that k − has a relatively small variation and the averages of k + and k − are the same (k), then even a highly variable distribution for k + , e.g., a bimodal distribution P(k + ) = [δ (k + ) + δ (k + − 2k)]/2 with a maximum fold change for k + , k +,max k +,min → ∞ would only lead to a modest fold change in k: f ≈ 3, which is far below what is needed to achieve V ≈ 2 within the 2-state model.

The model parameters are obtained with an optimization algorithm
Given that we can derive analytical expressions for the relevant quantities using the model equations, we implemented a gradient descent algorithm to find the rate parameters leading to the best agreement between the theory and the experimental results.
More precisely, we tried to find the values of k + (N), c i (i = 1, 2, 3), σ + and σ − , as defined in the previous sections, that would lead to values of τ , f + and V N = σ 2 τ / τ 2 N where ... N indicates an average over N, that are closest to what can be directly obtained from experiments. To do this, we define a vector parameter p, whose components are the parameters that we are trying to fit, and a loss function L defined as where α 1 , α 2 and α 3 are three constant weights. At every step the parameter vector is updated according to where η is a constant learning rate. Here, with exp and mod we indicate the mean values of each quantity as measured from the experimental data or as obtained from the model equations respectively.
The fit of our model to the experimental data shown in Fig. 3 in the main text led to the parameters c 1 = 0.30 ± 0.11, c 2 = 0.12 ± 0.14, c 3 = 0.06 ± 0.16, σ + = 0.19 ± 0.1, σ − ≤ 0.0005. The errors correspond to the fluctuation in the value of the parameter that would lead to a 50% increase in the loss function.
The parameters may also be obtained or further refined by using the maximal likelihood estimation (MLE), which fits the full distribution that includes in principle all the higher order moments of the distribution. Unless the model is inconsistent with the data, the "adjustment" of these parameters is typically small and they may not lead to a better fit of the physically meaningful properties of the system such as τ or f + .
We performed MLE for our problem, and indeed it does not change the results significantly. In particular, we performed the optimization again using a MLE loss function of the following form: where x are the experimental data, p the fitting parameter, P exp and P th the dwell time distributions for a given number of stator units, D is a measure of the distance between the two functions and dt the bin size we chose to evaluate the distributions. We used 8/12 the mean squared difference here for D. Other distance measures such as the log likelihood function (or the Kullback-Leibler divergence) lead to poor fitting of the distribution and τ because these distance measures tend to overemphasize the outliers in the dwell time distribution at large τ where the noise in the measured distribution is high.
We used the gradient descent method to minimize L MLE given above and the parameters obtained are only slightly different from those we reported in the article, namely c 1 = 0.39 (0.3 ± 0.1 in the article), c 2 = 0.14 (0.12 ± 0.14 in the article), c 3 = 0.012 (0.06 ± 0.16 in the article), σ + = 0.09 (0.19 ± 0.1 in the article). Although the predicted values for τ, f + , and V from these adjusted parameters agree with experiments within the experimental error bars, the agreement is not improved in comparison with those shown in Fig. 3 in the main text with parameters determined by minimizing the moment-based loss function of Eq. 28. Since the adjustment in parameters was small and the fitting to τ and f + is not improved, we did not use the MLE approach in our study. Supplementary Figure 4. a. Histograms of dwell times ending with binding (above) and unbinding (below) of a new stator unit, for three different values of N. b. Conditional dwell times τ + (yellow circle) and τ − (blue triangle) from the electrorotation experiments vs the number of stator units bound to the rotor. Data presented are the means and the error bars represent standard deviation. The sample size for determining τ + for 0, 1, 2,... 8 stator units is 74, 51, 54, 60, 54, 67, 52, 54, 37, respectively. The sample size for determining τ − for 1, 2,... 8 stator units is 12, 21, 10, 23, 14, 28, 29, 37, respectively.